function qr_iteration_test()
% QR_ITERATION_TEST QR迭代法专项测试
% 
% 详细测试QR迭代法的各种变形和优化技术

fprintf('=== QR迭代法专项测试 ===\n\n');

%% 测试1：基本QR迭代
fprintf('测试1：基本QR迭代\n');
fprintf('目标：理解基本QR迭代的收敛行为\n\n');

% 测试矩阵
A1 = [4, 1, 0; 1, 3, 1; 0, 1, 2];
fprintf('测试矩阵 A1 =\n');
disp(A1);

exact_eigenvals = sort(eig(A1), 'descend');
fprintf('精确特征值: [%.6f, %.6f, %.6f]\n', exact_eigenvals);

% 基本QR迭代
A_qr = A1;
fprintf('\n基本QR迭代过程:\n');
fprintf('迭代    对角元素                    下三角范数    收敛率\n');
fprintf('----    --------                    ----------    ------\n');

prev_off_diag = inf;
for k = 1:50
    [Q, R] = qr(A_qr);
    A_qr = R * Q;
    
    % 计算收敛指标
    off_diag_norm = norm(tril(A_qr, -1), 'fro');
    conv_rate = off_diag_norm / prev_off_diag;
    
    if k <= 15 || mod(k, 10) == 0 || off_diag_norm < 1e-12
        diag_elements = diag(A_qr);
        fprintf('%4d    [%6.4f, %6.4f, %6.4f]    %10.2e    %6.3f\n', ...
                k, diag_elements, off_diag_norm, conv_rate);
    end
    
    if off_diag_norm < 1e-12
        fprintf('基本QR迭代收敛于第 %d 次迭代\n', k);
        break;
    end
    
    prev_off_diag = off_diag_norm;
end

final_eigenvals = sort(diag(A_qr), 'descend');
fprintf('QR迭代结果: [%.6f, %.6f, %.6f]\n', final_eigenvals);
fprintf('误差: [%.2e, %.2e, %.2e]\n', abs(final_eigenvals - exact_eigenvals));

%% 测试2：带位移的QR迭代
fprintf('\n测试2：带位移的QR迭代\n');
fprintf('目标：比较不同位移策略的效果\n\n');

% 测试不同的位移策略
shift_strategies = {'无位移', 'Rayleigh商位移', 'Wilkinson位移'};
A_matrices = {A1, A1, A1};

for strategy_idx = 1:length(shift_strategies)
    fprintf('2.%d %s:\n', strategy_idx, shift_strategies{strategy_idx});
    
    A_shift = A_matrices{strategy_idx};
    iteration_count = 0;
    
    for k = 1:30
        n = size(A_shift, 1);
        
        % 选择位移
        switch strategy_idx
            case 1  % 无位移
                shift = 0;
            case 2  % Rayleigh商位移
                shift = A_shift(n, n);
            case 3  % Wilkinson位移
                % 计算右下角2×2子矩阵的特征值，选择更接近A(n,n)的
                if n >= 2
                    sub_matrix = A_shift(n-1:n, n-1:n);
                    sub_eigenvals = eig(sub_matrix);
                    [~, idx] = min(abs(sub_eigenvals - A_shift(n, n)));
                    shift = sub_eigenvals(idx);
                else
                    shift = A_shift(n, n);
                end
        end
        
        % QR分解和更新
        [Q, R] = qr(A_shift - shift * eye(n));
        A_shift = R * Q + shift * eye(n);
        
        off_diag = norm(tril(A_shift, -1), 'fro');
        iteration_count = k;
        
        if k <= 8 || off_diag < 1e-12
            fprintf('  迭代 %2d: 位移 = %8.4f, 收敛指标 = %.2e\n', ...
                    k, shift, off_diag);
        end
        
        if off_diag < 1e-12
            break;
        end
    end
    
    fprintf('  收敛迭代数: %d\n', iteration_count);
    final_vals = sort(diag(A_shift), 'descend');
    fprintf('  最终特征值: [%.6f, %.6f, %.6f]\n\n', final_vals);
end

%% 测试3：Hessenberg预处理
fprintf('测试3：Hessenberg预处理\n');
fprintf('目标：测试Hessenberg化对QR迭代的加速效果\n\n');

% 构造较大的测试矩阵
n3 = 6;
A3 = randn(n3, n3);
A3 = (A3 + A3') / 2;  % 对称化以确保实特征值

fprintf('测试矩阵规模: %d×%d\n', n3);
exact_eigenvals3 = sort(eig(A3), 'descend');

% 直接QR迭代
fprintf('\n3.1 直接QR迭代:\n');
A_direct = A3;
tic;
for k = 1:200
    [Q, R] = qr(A_direct);
    A_direct = R * Q;
    
    if norm(tril(A_direct, -1), 'fro') < 1e-10
        break;
    end
end
time_direct = toc;
fprintf('直接QR迭代: %d次迭代, %.4f秒\n', k, time_direct);

% Hessenberg预处理 + QR迭代
fprintf('\n3.2 Hessenberg预处理 + QR迭代:\n');
tic;
[P, A_hess] = hess(A3);  % Hessenberg化
hess_time = toc;

fprintf('Hessenberg化时间: %.4f秒\n', hess_time);
fprintf('Hessenberg矩阵:\n');
disp(A_hess);

tic;
for k = 1:200
    [Q, R] = qr(A_hess);
    A_hess = R * Q;
    
    if norm(tril(A_hess, -2), 'fro') < 1e-10  % 注意：Hessenberg矩阵下三角从-2开始
        break;
    end
end
time_hess_qr = toc;

total_time_hess = hess_time + time_hess_qr;
fprintf('Hessenberg QR迭代: %d次迭代, %.4f秒\n', k, time_hess_qr);
fprintf('总时间: %.4f秒\n', total_time_hess);
fprintf('加速比: %.2f\n', time_direct / total_time_hess);

%% 测试4：双重位移QR算法
fprintf('\n测试4：双重位移QR算法\n');
fprintf('目标：处理复特征值的情况\n\n');

% 构造有复特征值的矩阵
theta = pi/6;
A4 = [cos(theta), -sin(theta), 0.1; sin(theta), cos(theta), 0.2; 0, 0, 0.5];
fprintf('有复特征值的矩阵 A4:\n');
disp(A4);

exact_eigenvals4 = eig(A4);
fprintf('精确特征值:\n');
for i = 1:length(exact_eigenvals4)
    if imag(exact_eigenvals4(i)) == 0
        fprintf('  λ%d = %.6f\n', i, real(exact_eigenvals4(i)));
    else
        fprintf('  λ%d = %.6f %+.6fi\n', i, real(exact_eigenvals4(i)), imag(exact_eigenvals4(i)));
    end
end

% 实Schur形式
fprintf('\n使用实Schur分解:\n');
[T, Z] = schur(A4, 'real');
fprintf('实Schur形式 T:\n');
disp(T);

% 从Schur形式提取特征值
schur_eigenvals = [];
i = 1;
while i <= size(T, 1)
    if i == size(T, 1) || abs(T(i+1, i)) < 1e-12
        % 实特征值
        schur_eigenvals = [schur_eigenvals; T(i, i)];
        i = i + 1;
    else
        % 复共轭特征值对
        sub_block = T(i:i+1, i:i+1);
        complex_eigenvals = eig(sub_block);
        schur_eigenvals = [schur_eigenvals; complex_eigenvals];
        i = i + 2;
    end
end

fprintf('从Schur形式提取的特征值:\n');
for i = 1:length(schur_eigenvals)
    if imag(schur_eigenvals(i)) == 0
        fprintf('  λ%d = %.6f\n', i, real(schur_eigenvals(i)));
    else
        fprintf('  λ%d = %.6f %+.6fi\n', i, real(schur_eigenvals(i)), imag(schur_eigenvals(i)));
    end
end

%% 测试5：收敛性分析
fprintf('\n测试5：收敛性分析\n');
fprintf('目标：分析不同矩阵类型的收敛行为\n\n');

matrix_types = {
    '对角占优', diag([5, 3, 1]) + 0.1*randn(3,3);
    '近似对角', diag([4, 2, 1]) + 0.5*randn(3,3);
    '一般矩阵', randn(3,3);
    '病态矩阵', [1, 1, 1; 0, 1e-6, 0; 0, 0, 1e-12]
};

fprintf('矩阵类型      条件数      收敛迭代数    最终精度\n');
fprintf('--------      ------      ----------    --------\n');

for i = 1:size(matrix_types, 1)
    type_name = matrix_types{i, 1};
    A_test = matrix_types{i, 2};
    
    % 确保矩阵是实对称的以避免复特征值
    A_test = (A_test + A_test') / 2;
    
    cond_num = cond(A_test);
    
    % QR迭代
    A_conv = A_test;
    for k = 1:500
        [Q, R] = qr(A_conv);
        A_conv = R * Q;
        
        off_diag = norm(tril(A_conv, -1), 'fro');
        if off_diag < 1e-10
            break;
        end
    end
    
    final_precision = norm(tril(A_conv, -1), 'fro');
    
    fprintf('%-12s  %8.2e    %10d    %8.2e\n', ...
            type_name, cond_num, k, final_precision);
end

%% 测试6：大规模性能测试
fprintf('\n测试6：大规模性能测试\n');
fprintf('目标：测试QR迭代在大规模矩阵上的性能\n\n');

sizes = [20, 50, 100];
fprintf('矩阵规模    QR迭代时间    Hessenberg+QR    MATLAB eig\n');
fprintf('--------    ----------    -------------    ----------\n');

for n = sizes
    % 生成随机对称矩阵
    A_large = randn(n, n);
    A_large = (A_large + A_large') / 2;
    
    % 直接QR迭代（限制迭代次数以避免过长时间）
    tic;
    A_qr_large = A_large;
    max_iter = min(100, n*2);  % 限制迭代次数
    for k = 1:max_iter
        [Q, R] = qr(A_qr_large);
        A_qr_large = R * Q;
        
        if norm(tril(A_qr_large, -1), 'fro') < 1e-8
            break;
        end
    end
    time_qr = toc;
    
    % Hessenberg + QR
    tic;
    [P, A_hess_large] = hess(A_large);
    for k = 1:max_iter
        [Q, R] = qr(A_hess_large);
        A_hess_large = R * Q;
        
        if norm(tril(A_hess_large, -2), 'fro') < 1e-8
            break;
        end
    end
    time_hess = toc;
    
    % MATLAB内置eig
    tic;
    eigenvals_matlab = eig(A_large);
    time_matlab = toc;
    
    fprintf('%6d      %10.4f    %13.4f    %10.4f\n', ...
            n, time_qr, time_hess, time_matlab);
end

%% 测试7：特殊结构矩阵
fprintf('\n测试7：特殊结构矩阵\n');
fprintf('目标：测试QR迭代在特殊结构矩阵上的表现\n\n');

% 三对角矩阵
fprintf('7.1 三对角矩阵:\n');
n7 = 5;
A_tridiag = diag(2*ones(n7,1)) + diag(ones(n7-1,1), 1) + diag(ones(n7-1,1), -1);
fprintf('三对角矩阵:\n');
disp(A_tridiag);

A_tri_qr = A_tridiag;
for k = 1:50
    [Q, R] = qr(A_tri_qr);
    A_tri_qr = R * Q;
    
    if norm(tril(A_tri_qr, -1), 'fro') < 1e-12
        break;
    end
end

fprintf('QR迭代收敛: %d次迭代\n', k);
tri_eigenvals = sort(diag(A_tri_qr), 'descend');
exact_tri_eigenvals = sort(eig(A_tridiag), 'descend');
fprintf('特征值误差: %.2e\n', norm(tri_eigenvals - exact_tri_eigenvals));

% 上Hessenberg矩阵
fprintf('\n7.2 上Hessenberg矩阵:\n');
A_hess_test = triu(randn(4,4), -1);  % 上Hessenberg
fprintf('上Hessenberg矩阵:\n');
disp(A_hess_test);

A_hess_qr = A_hess_test;
for k = 1:30
    [Q, R] = qr(A_hess_qr);
    A_hess_qr = R * Q;
    
    if norm(tril(A_hess_qr, -2), 'fro') < 1e-12
        break;
    end
end

fprintf('Hessenberg QR收敛: %d次迭代\n', k);

fprintf('\n=== QR迭代法测试总结 ===\n');
fprintf('1. 基本QR迭代收敛较慢，但最终能达到高精度\n');
fprintf('2. 位移策略显著加速收敛，Wilkinson位移效果最好\n');
fprintf('3. Hessenberg预处理大幅提高计算效率\n');
fprintf('4. 实Schur分解可以处理复特征值情况\n');
fprintf('5. 矩阵条件数影响收敛速度和精度\n');
fprintf('6. 对于大规模问题，预处理是必要的\n');
fprintf('7. 特殊结构矩阵可以利用结构加速计算\n\n');

end